*This program runs a probit and calculates the associated marginal effects
*for stocks, own business and home in ELSA
*It's used on the paper on wealth counterfactuals 
*for the final verification for the Review of Economics and Statistics
*June 2011 - Dimitris Christelis

version 9.2
cap clear all
macro drop _all
clear
set mem 500m
set more off

*Paths
do "c:/foldp/foldp.do"
global projf "${foldp}/Wealth Counterfactuals/Verification"

*Tolerance
sca tolr=1e-6

*Weights
global wgt wgth

*ATTENTION: The base country (Germany in our case) has to be put first in the sequence
global eureg1 uk
global eureg2 England
*Country numbers: they have to follow exactly the same sequence as the country names
global eureg3 12   
local nas: word count $eureg1

global assetl1 rfin home ownb secdebt liab tdebt
global assetl2 hrfin hhome hownb hsecdebt hliab htdebt
local nar: word count $assetl1


*Attention: the asset loop has to be before the country loop
forvalues ll=1/`nar' {

   *Globals of asset
   global as: word `ll' of $assetl1
   global hhas: word `ll' of $assetl2

   forvalues kk=1/`nas' {

   *Globals of country
   global re: word `kk' of $eureg1
   global fnre: word `kk' of $eureg2
   global cn: word `kk' of $eureg3


   *Logging
   cap log close
   log using "${projf}/Programs/ELSA/${re}_${as}.log", replace

   global pdate "$S_DATE" 
   global ptime "$S_TIME"
   display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"

   *Opening HRS data to generate quartiles
   use "${projf}/Data/HRS/HRS_DG.dta", clear

   *Checking the sample weights
   assert abs(wgth-wgtach)<tolr
   drop wgtach

   *Keeping only heads
   keep if head==1 & wgth>0 & wgth!=.

   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov /*+ hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Creating wealth quartiles 
   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_us = r(p25)
   sca p50_w_us = r(p50)
   sca p75_w_us = r(p75)

   *Creating income quartiles
   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_us = r(p25)
   sca p50_i_us = r(p50)
   sca p75_i_us = r(p75)

   *Reading the SHARE data
   use "${projf}/Data/SHARE/SHARE_DG.dta", clear

   *Keeping only the country of interest
   keep if country==3 & ${wgt}>0 & ${wgt}!=.

   *Keeping only heads
   keep if head==1

   *Total debt
   qui gen double htdebtv = hmortv+hliabv 
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt
   qui ge double hsecdebtv = hmortv 
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_de = r(p25)
   sca p50_w_de = r(p50)
   sca p75_w_de = r(p75)

   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_de = r(p25)
   sca p50_i_de = r(p50)
   sca p75_i_de = r(p75)

   *Reading the parameters, variable lists etc.
   do "${projf}/Programs/Parameters.do"

   use "${projf}/Data/ELSA/ELSA_DG.dta", clear

   *Keeping only the country of interest
   keep if country==${cn} & ${wgt}>0 & ${wgt}!=.

   *Keeping only heads
   keep if head==1

   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov /*+ hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

  
   *Rescaling age
   qui replace age=age/sc_age
   qui gen double agesq=age^2

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   qui gen double nw=hnetwv-${hhas}v
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Bad health dummy
   qui gen byte badhs=(srhealtha>=4)

   *College dummy
   qui gen byte college=(edu==3)
   
   *High school dummy
   qui gen byte highs=(edu==2)

   *Couple dummy
   cap drop couple
   qui gen byte couple=(mstatus==1)
   
   *Widow dummy
   qui gen byte widow=(mstatus==3)
   
   *Never married
   qui gen byte nevmar=(mstatus==4)
      
   *Creating wealth quartiles for financial and real assets
   qui gen double wqnw1_qus=(nw<=scalar(p25_w_us))
   qui gen double wqnw2_qus=(nw>scalar(p25_w_us) & nw<=scalar(p50_w_us))
   qui gen double wqnw3_qus=(nw>scalar(p50_w_us) & nw<=scalar(p75_w_us))
   qui gen double wqnw4_qus=(nw>scalar(p75_w_us))
   qui gen double wqnw1_qde=(nw<=scalar(p25_w_de))
   qui gen double wqnw2_qde=(nw>scalar(p25_w_de) & nw<=scalar(p50_w_de))
   qui gen double wqnw3_qde=(nw>scalar(p50_w_de) & nw<=scalar(p75_w_de))
   qui gen double wqnw4_qde=(nw>scalar(p75_w_de))

   qui gen double wqinc1_qus=(inc<=scalar(p25_i_us))
   qui gen double wqinc2_qus=(inc>scalar(p25_i_us) & inc<=scalar(p50_i_us))
   qui gen double wqinc3_qus=(inc>scalar(p50_i_us) & inc<=scalar(p75_i_us))
   qui gen double wqinc4_qus=(inc>scalar(p75_i_us))
   qui gen double wqinc1_qde=(inc<=scalar(p25_i_de))
   qui gen double wqinc2_qde=(inc>scalar(p25_i_de) & inc<=scalar(p50_i_de))
   qui gen double wqinc3_qde=(inc>scalar(p50_i_de) & inc<=scalar(p75_i_de))
   qui gen double wqinc4_qde=(inc>scalar(p75_i_de))

   *Saving the dataset
   qui save "${projf}/Results/ELSA/${as}/temp.dta", replace

   *Running 2 regressions, using both definitions of quartiles
   foreach cb in qus qde {

      use "${projf}/Results/ELSA/${as}/temp.dta", clear
      probit ${hhas}o ${basiccovs_`cb'}, r

      *Coefficient matrix
      mat olco=e(b)
  
      *Augmenting the coefficient matrix by an empty column, the number of obs and the log likelihood
      mat blank = .
      mat colnames blank=blank
  
      mat olN = e(N)
      mat colnames olN=nobs
  
      mat ill = e(ll)
      mat colnames ill=loglik

      mat olco=[olco,blank,olN,ill]
       
      *Var-Covar matrix
      mat olvar=e(V)

      *Saving matrix of regression results as a text file
  
      *Number of variables
      local nc=colsof(olvar)
      global nc=`nc'
      matrix olout=J(`nc'+3,6,.)

      *Putting in the regression coefficients 
      forvalues k=1/`nc' {
         matrix olout[`k',1]=olco[1,`k']
         matrix olout[`k',2]=sqrt(olvar[`k',`k'])
         matrix olout[`k',3]=olco[1,`k']/sqrt(olvar[`k',`k'])
         matrix olout[`k',4]=2*ttail(e(N) - `nc',abs(olout[`k',3]))
         matrix olout[`k',5] = olout[`k',1] - invttail(e(N) - `nc',sigl/2)*olout[`k',2]
         matrix olout[`k',6] = olout[`k',1] + invttail(e(N) - `nc',sigl/2)*olout[`k',2]
      }
  
      *Adding number of observations used in first stage probit
      matrix olout [`nc'+2,1]=e(N)
      *Adding log likelihood
      matrix olout [`nc'+3,1]=e(ll)

      *Labels of matrix of regression output
      *Column labels
      local colmat coeff se tstat pvalue ci_low ci_high
      matrix colnames olout = `colmat' 
      *Row labels
      matrix rownames olout = ${basiccovs_`cb'} _cons blank nobs loglik
      *Displaying and saving the matrix of results
      mat list olout,format(%15.5f)

      mat2txt, mat(olout) saving("${projf}/Results/ELSA/${as}/probit_est_${re}_${as}_`cb'.txt") ///
         replace format(%15.5f) title("Regression for ${as} in ${fnre}, using `cb' quartiles")
  
      preserve
      *Saving matrices of coefficients and correlations and their covariance matrices
      *for each implicate
      drop _all
      svmat double olco, names(eqcol)
      qui save "${projf}/Results/ELSA/${as}/probit_coeff_${re}_${as}_`cb'.dta", replace
      drop _all
      svmat double olvar, names(eqcol)
      qui save "${projf}/Results/ELSA/${as}/probit_var_${re}_${as}_`cb'.dta", replace
      restore  

      *Bootstrap program

      *Creating a matrix with number of columsn equal to 1 + # bootstrap runs, storing
      *the main sample estimates in the first line, together with the proportion of
      *hhds owning the asset
      qui su ${hhas}o if head==1 [aw=$wgt], meanonly
      sca ${as}p=r(mean)
      mat olco=olco[1,1..${nc}]
      mat ${as}coeffall=(olco,${as}p)
      local nch=$nc+1
      local bn ""
      forvalues lll=1/`nch' {
        local bn `bn' b`lll'
      }
      mat colnames ${as}coeffall = `bn'   
      qui save "${projf}/Results/ELSA/${as}/temp2.dta", replace
  
      *Subsequent runs, beginning of the bootstrap loop
      forvalues i=1/$finsim {

         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
           
         use "${projf}/Results/ELSA/${as}/temp2.dta", clear
         *Fixing the seed before each resampling
         local sd=1234+30*`i'+`kk'
         set seed `sd'
         bsample
         
         *Probit regression in the bootstrapped sample
         if `i'==1 {
            probit ${hhas}o ${basiccovs_`cb'}, r
         }
         if `i'>1 {
            qui probit ${hhas}o ${basiccovs_`cb'}, r
         }

         mat a=e(V)
         local nc`i'=colsof(a)
               
         *If no variables are dropped from the probit then we proceed with updating the matrix
         if `nc`i''==$nc {
            qui su ${hhas}o if head==1 [aw=$wgt], meanonly
            sca ${as}p=r(mean)
            mat olco=e(b)
            mat ${as}coeff`sr'=(olco,${as}p)
            mat ${as}coeffall=(${as}coeffall \ ${as}coeff`sr')

            forvalues lll=1/10 {
               if `sr'==$nsim*(`lll'/10) {
                  noisily di in yellow `sr' in yellow "*" _continue 
               }
            }

            mat drop ${as}coeff`sr'
         }
         
         *If some variables are dropped from the probit in the bootstrapped samples 
         *we put back the index to its previous value
         if `nc`i''~=$nc {
            local sr=`sr'-1
         }
               
         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         if `sr'==$nsim {
            di in yellow "End of ${nsim} valid boostrap runs, total runs equal to `i'"
            continue, break
         }
         
      }  /* End of loop over bootstrap runs */

      *Saving the matrix of bootstrap coefficients and sample proportions
      drop _all
      svmat double ${as}coeffall, names(eqcol)
      qui save "${projf}/Results/ELSA/${as}/boot_drawcoeff_${re}_${as}_`cb'.dta", replace


***********************************************************
*Calculation of marginal effects and associated magnitudes*
***********************************************************

      *Reading the matrices saved from the estimation part, to retrieve the # of obs
      use "${projf}/Results/ELSA/${as}/probit_coeff_${re}_${as}_`cb'.dta", clear
      mkmat _all , mat(a) 
      local nc=colsof(a)
      global nc=`nc'-3
      sca nobs_`cb'=a[1,${nc}+2]


      *Reading the dataset
      use "${projf}/Results/ELSA/${as}/temp.dta", clear

      *Incrementing the continuous variables by the appropriate magnitudes
      *to compute marginal effects and elasticities
      *Variables that are not used in elasticity form
      foreach y in $conlist_meff {
         qui gen double `y'_pl=`y'+(ch_`y'/sc_`y')
      }   

      merge using "${projf}/Results/ELSA/${as}/boot_drawcoeff_${re}_${as}_`cb'.dta"   
      drop _merge   

      *Putting the values of the drawn coefficients equal to the dummy value from missing
      *in rows lower than the number of bootstraps, so that there are no missing values
      local nch = $nc + 1
      forvalues h=1/`nch' {   
         qui replace _b`h'=$dumval if _n>${nsim}+1  
      }      
   
      *Renaming the coefficients  
      local ncl = $nc - 1
      forvalues h=1/`ncl' {   
         local g: word `h' of ${basiccovs_`cb'}
         rename _b`h' b_`g'_`cb'   
      }   

      *Renaming the coefficient of the constant 
      *set trace on
      rename _b${nc} b_cons_`cb'
   
      local nch = $nc + 1
      rename _b`nch' pr_${re}_${as}_`cb'_own
   
      *Merging with the draws from the US
      if "`cb'"=="qus" {
         merge using "${projf}/Results/HRS/Country/${as}/boot_drawcoeff_usall_${as}.dta"   
      	 drop _merge   

      	 *Making sure that the number of coefficients in the ELSA regression is the 
      	 *same as the one in the US regression
      	 preserve
      	 use "${projf}/Results/HRS/Country/${as}/probit_var_usall_${as}.dta", clear
      	 mkmat _all , mat(olvar) 
      	 local chnc=colsof(olvar)
      	 assert `chnc'==${nc}
      	 restore

      	 *Global of number of variables
      	 global nc=colsof(olvar)


      	 *Putting the values of the drawn coefficients equal to the dummy value from missing
      	 *in rows lower than the number of bootstraps, so that there are no missing values
      	 local nch = $nc + 1
      	 forvalues h=1/`nch' {   
      	    qui replace _b`h'=$dumval if _n>${nsim}+1  
      	 }      
   
      	 *Renaming the coefficients  
      	 local ncl = $nc - 1
      	 forvalues h=1/`ncl' {   
      	    local g: word `h' of ${basiccovs_qus}
      	    rename _b`h' b_usall_`g'   
      	 }   

      	 *Renaming the coefficient of the constant 
      	 *set trace on
      	 rename _b${nc} b_usall_cons
   
      	 local nch = $nc + 1
      	 rename _b`nch' pr_usall_${as}_own
      }
      
      *Merging with the coefficients and the probabilities from Germany
      if "`cb'"=="qde" {
         if "${re}"~="de" {

      	    *Coefficients
      	    merge using "${projf}/Results/SHARE/de/${as}/boot_drawcoeff_de_${as}_qde.dta"   
      	    drop _merge   

      	    *Putting the values of the drawn coefficients equal to the dummy value from missing
      	    *in rows lower than the number of bootstraps, so that there are no missing values
      	    local nch = $nc + 1
      	    forvalues h=1/`nch' {   
      	       qui replace _b`h'=$dumval if _n>${nsim}+1  
      	    }      
   
      	    *Renaming the coefficients  
      	    local ncl = $nc - 1
      	    forvalues h=1/`ncl' {   
      	       local g: word `h' of ${basiccovs_qde}
      	       rename _b`h' b_de_`g'_qde   
      	    }   

      	    *Renaming the coefficient of the constant 
      	    *set trace on
      	    rename _b${nc} b_de_cons_qde
   
      	    local nch = $nc + 1
      	    rename _b`nch' pr_de_${as}_qde_own
      	 }         
      }
   
      *Loop over bootstrap coefficients
      forvalues i=1/$finsim { 

         *Defining an index of valid bootstrap runs that will determine whether the 
         *loop finishes or not
         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
         *If there were any condition on the coefficients that was violated by the current
         *draw then the index of valid bootstrap runs would be decreased by one at this point
         *in the program. For the pooled probit there is no such condition and thus the
         *index of valid bootstrap runs (sr) is equal to the on of total bootstrap runs (i)
   
         *Creating the linear prediction for the current draw of the coefficients
         qui gen double xb0_`cb' = b_cons_`cb'[`i']
         local lll ${basiccovs_`cb'}
         foreach y in `lll' {
            qui replace xb0_`cb' = xb0_`cb' + b_`y'_`cb'[`i']*`y' 
         }
      
         *Creating a counterfactual prediction using the coefficients from the US
         if "`cb'"=="qus" {
            qui gen double xb0_usall = b_usall_cons[`i']
            foreach y in $basiccovs_qus {
               qui replace xb0_usall = xb0_usall + b_usall_`y'[`i']*`y' 
            }
	 }

         *Creating a counterfactual prediction using the coefficients from Germany
         if "`cb'"=="qde" {
            if "${re}"~="de" {         
               qui gen double xb0_de_qde = b_de_cons_qde[`i']
               foreach y in $basiccovs_qde {
                  qui replace xb0_de_qde = xb0_de_qde + b_de_`y'_qde[`i']*`y' 
               }
            }            
         }
      
         *******************************************************************
         *Calculating the average of probabilities across individual obs
         if "`cb'"=="qus" {   
            qui gen double dp_usall=normal(xb0_usall)
            qui su dp_usall [aw=$wgt], meanonly
            sca r0m=r(mean)
         }
         
         *Probability using Germany coefficients
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               qui gen double dp_de_qde=normal(xb0_de_qde)
               qui su dp_de_qde [aw=$wgt], meanonly
               sca r2m_qde=r(mean)
            }
         }
*set trace on
         if `i'==1 {
             
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients
               qui gen double pr_${re}_${as}_qus_cus=r0m
             
               *Difference between US ownership proportion and
               *own ownership proportion
               qui gen double dpus_${re}_${as}_qus_tot=pr_usall_${as}_own - ///
                  pr_${re}_${as}_qus_own
             
               *Difference between US and own probability calculated using
               *US coefficients (characteristics effect)
               qui gen double dpus_${re}_${as}_qus_char=pr_usall_${as}_own[1] - ///   
                  pr_${re}_${as}_qus_cus[1]
             
               *Difference between own probability calculated using
               *US coefficients and own probability (coefficient effect)
               qui gen double dpus_${re}_${as}_qus_coeff=pr_${re}_${as}_qus_cus[1] - ///   
                  pr_${re}_${as}_qus_own[1]
            } 

            if "`cb'"=="qde" {
               if "${re}"~="de" {
            
             	  *Counterfactual probability using German coefficients
             	  qui gen double pr_${re}_${as}_qde_cde=r2m_qde
             
             	  *Difference between mean German predicted probability and
             	  *mean own predicted probability
             	  qui gen double dpde_${re}_${as}_qde_tot=pr_de_${as}_qde_own - ///
             	     pr_${re}_${as}_qde_own
             
             	  *Difference between German and own probability calculated using
             	  *German coefficients (characteristics effect)
             	  qui gen double dpde_${re}_${as}_qde_char = pr_de_${as}_qde_own[1] - ///
             	     pr_${re}_${as}_qde_cde[1]  
             
             	  *Difference between own probability calculated using
             	  *German coefficients and own probability (coefficient effect)
             	  qui gen double dpde_${re}_${as}_qde_coeff = pr_${re}_${as}_qde_cde[1] - ///
             	    pr_${re}_${as}_qde_own[1]
               }
            }
         }
         
         *Putting the probability at the line indexed by
         *the valid bootstrap run (sr), and not the general
         *run (i)
         
         if `i'>1 {
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients   
               qui replace pr_${re}_${as}_qus_cus=r0m in `sr'   
               assert pr_${re}_${as}_qus_cus<$dumval if _n==`sr'   

               *Difference between US and own probability calculated using   
               *US coefficients (characteristics effect)   
               qui replace dpus_${re}_${as}_qus_char=pr_usall_${as}_own[`sr'] - ///   
                  pr_${re}_${as}_qus_cus[`sr'] in `sr'   

               *Difference between own probability calculated using   
               *US coefficients and own probability (coefficient effect)   
               qui replace dpus_${re}_${as}_qus_coeff=pr_${re}_${as}_qus_cus[`sr'] - ///   
                  pr_${re}_${as}_qus_own[`sr'] in `sr'   
            }

            if "`cb'"=="qde" {
               if "${re}"~="de" {
                  *Counterfactual probability using German coefficients
               	  qui replace pr_${re}_${as}_qde_cde=r2m_qde in `sr'
               	  assert pr_${re}_${as}_qde_cde<$dumval if _n==`sr'
         
               	  *Difference between German and own probability calculated using
               	  *German coefficients (characteristics effect)
               	  qui replace dpde_${re}_${as}_qde_char=pr_de_${as}_qde_own[`sr'] - ///
               	     pr_${re}_${as}_qde_cde[`sr'] in `sr'

               	  *Difference between own probability calculated using
               	  *German coefficients and own probability (coefficient effect)
               	  qui replace dpde_${re}_${as}_qde_coeff=pr_${re}_${as}_qde_cde[`sr'] - ///
               	     pr_${re}_${as}_qde_own[`sr'] in `sr'
               } 
            }
         }
         
         if "`cb'"=="qus" {
            drop dp_usall xb0_usall
         }
         
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               drop dp_de_qde xb0_de_qde
            }   
         }   

         
         ******************************************************************************************        	 
         ************* 	
         *0/1 Dummies* 	
         ************* 	
    	
         foreach y in $dumlist_meff /* $dumlist */ {           	

            *Putting all dummies to zero          	
            *Removing the effect of the dummy.         	
            qui gen double xb_c1_`cb' = xb0_`cb' - b_`y'_`cb'[`i']*`y'    	        	
            assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
            *Putting back the dummy equal to 1 	
            qui gen double xb_c2_`cb' = xb_c1_`cb' + b_`y'_`cb'[`i']    	        	
            assert xb_c2_`cb'~=. if xb0_`cb'~=.   	        	
    	
            *Loop over the marginal effect and percentage variation 	
            forvalues tt=1/$wcp { 	
               local wpf: word `tt' of ${prfx} 	
     	
      	       if `tt'==1 {    	
                  qui gen double pdif_`y'_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
               }    	
    	       if `tt'==2 {             	
	          qui ge double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	             normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
	       }    	
	       assert pdif_`y'_`cb'~=. if xb0_`cb'~=.   	
               qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                    	

   	       if `i'==1 {                   	
	          qui ge double `wpf'_`y'_`cb' = $dumval                   	
	       }                   	
               *Putting the effect at the appropriate line in the file
               qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                      	
	       assert `wpf'_`y'_`cb'<$dumval if _n==`sr' 	
               drop pdif_`y'_`cb'
	     		         	
            }   /* End of loop over magnitudes */ 	
            drop xb_c1_`cb' xb_c2_`cb' 	
         }  /* End of loop over list of 0/1 dummies */        	
	 
         
         ******************************************************************************************        	 
         ************************ 	
         *Interdependent Dummies* 	
         ************************ 	
         foreach g in ${groupno_meff_`cb'} /* forvalues g=1/$ng */ {  

            *Putting all dummies in a group equal to zero            
            qui gen double xb_c1_`cb' = xb0_`cb'
            local m ${group`g'}
            foreach y in `m'  {            
               qui replace xb_c1_`cb' = xb_c1_`cb' - b_`y'_`cb'[`i']*`y'            
            }	
            
            local m ${group`g'}
            foreach y in `m'  {     

               *Putting back the dummy equal to 1
     	       qui gen double xb_c2_`cb'= xb_c1_`cb' + b_`y'_`cb'[`i'] 

               *Loop over the elasticity and the 4 intervals
               forvalues tt=1/$wcp {
                  local wpf: word `tt' of ${prfx}    
          
	          if `tt'==1 {   
                     qui gen double pdif_`y'_`cb' = normal(xb_c2_`cb') - normal(xb_c1_`cb')     
	          }   
    	          if `tt'==2 {   
	             qui gen double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	                normal(xb_c1_`cb'))/normal(xb_c1_`cb')                     
	          }   
  	          assert pdif_`y'_`cb'~=. if xb0_`cb'~=.

  	          qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                     

  	          if `i'==1 {                     
	             qui ge double `wpf'_`y'_`cb' = $dumval                     
	          }                     
                  qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                     
	          assert `wpf'_`y'_`cb'<$dumval if _n==`sr'   
                  drop pdif_`y'_`cb' 
  
              }  /* End of loop over magnitudes */
               drop xb_c2_`cb' 
            }  /* End of loop over variables of a given group */       
            drop xb_c1_`cb'
         }  /* End of loop over groups */       

         ******************************************************************************************        	 
         ********************** 	
         *Continuous variables* 	
         ********************** 	
    	
         foreach y in $conlist_meff /* $conlist */ {           	

            qui gen double xb_c1_`cb' = xb0_`cb' 
            assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
            *Incrementing by the already determined sum	
            qui gen double xb_c2_`cb' = xb_c1_`cb' + b_`y'_`cb'[`i']*(`y'_pl - `y')    	        	
            assert xb_c2~=. if xb0_`cb'~=.   	        	
    	
            *Loop over the marginal effect and percentage variation 	
            forvalues tt=1/$wcp { 	
               local wpf: word `tt' of ${prfx} 	
     	
      	       if `tt'==1 {    	
                    qui gen double pdif_`y'_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
               }    	
    	       if `tt'==2 {             	
	          qui ge double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	             normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
	       }    	
	       assert pdif_`y'_`cb'~=. if xb0_`cb'~=.   	
               qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                    	
   	       if `i'==1 {                   	
	          qui ge double `wpf'_`y'_`cb' = $dumval                   	
	       }                   	
               *Putting the effect at the appropriate line in the file
               qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                      	
	       assert `wpf'_`y'_`cb'<$dumval if _n==`sr' 	
  	       drop pdif_`y'_`cb'	         	

            }   /* End of loop over magnitudes */ 	
            drop xb_c1_`cb' xb_c2_`cb' 	
         }  /* End of loop over list of continuous variables */        	
         
         ******************************************************************************************        	 
         *****
         *Age* 	
         *****

         qui gen double xb_c1_`cb' = xb0_`cb' 
         assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
         *Incrementing by the already determined sum	
         qui gen double xb_c2_`cb' = xb_c1_`cb' + b_age_`cb'[`i']*(ch_age/sc_age) + ///
            b_agesq_`cb'[`i']*(2*age*(ch_age/sc_age) + ((ch_age/sc_age)^2))
         assert xb_c2_`cb'~=. if xb0_`cb'~=.   	        	
    	
         *Loop over the marginal effect and percentage variation 	
         forvalues tt=1/$wcp { 	
            local wpf: word `tt' of ${prfx} 	
     	
            if `tt'==1 {    	
                 qui gen double pdif_age_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
            }    	
            if `tt'==2 {             	
               qui ge double pdif_age_`cb' = (normal(xb_c2_`cb') - ///
                 normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
            }    	
            assert pdif_age_`cb'~=.    	
            qui sum pdif_age_`cb' [aw = $wgt], meanonly                    	
            if `i'==1 {                   	
               qui ge double `wpf'_age_`cb' = $dumval                   	
            }                   	
            *Putting the effect at the appropriate line in the file
            qui replace `wpf'_age_`cb' = r(mean) in `sr'                      	
	    assert `wpf'_age_`cb'<$dumval if _n==`sr' 	
            drop pdif_age_`cb'	     		         	

         }   /* End of loop over magnitudes */ 	
         drop xb_c1_`cb' xb_c2_`cb' 	

         **********************************************************************************	 
         *Dropping linear prediction based on current draw
         drop xb0_`cb' 
      
         forvalues lll=1/10 {
            if `sr'==$nsim*(`lll'/10) {
               di in yellow "Iteration `sr', ${as}"
            }
         }

         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         local ncpp=$nsim+1
         if `sr'==`ncpp' {
            di in yellow "End of ${nsim} valid bootstrap runs, total runs equal to `i'"
            continue, break
         }
      
      }  /* End of loop over bootstrap runs */

      *Dropping obs beyond the line of the last valid bootstrap run
      drop if _n>$nsim+1
   
      *Global of marginal effects, putting first the probabilities of interest
      if "`cb'"=="qus" {
         local me pr_${re}_${as}_qus_own pr_${re}_${as}_qus_cus dpus_${re}_${as}_qus_tot ///
            dpus_${re}_${as}_qus_coeff dpus_${re}_${as}_qus_char
      }
      if "`cb'"=="qde" {
         if "${re}"=="de" {
            local me pr_${re}_${as}_qde_own
         } 
         if "${re}"~="de" {
            local me pr_${re}_${as}_qde_own pr_${re}_${as}_qde_cde dpde_${re}_${as}_qde_tot ///
               dpde_${re}_${as}_qde_coeff dpde_${re}_${as}_qde_char
         } 
      }
      
      foreach y in $dumlist_meff ${grouplist_meff_`cb'} $conlist_meff age {
         if "`y'"~="agesq" {
            local me `me' me_`y'
         }
      }
      
      *Keeping only the marginal effects and percentage variation
      keep country `me' `perc' 

      *Calculating the variance of the generated average probabilities and marginal effects
      *One has to start from the 2nd line, since the first line gives the magnitudes
      *evaluated at the estimated coefficients
      local ncpp=${nsim}+1
      foreach vr in `me' `perc' {
         qui su `vr' in 2/`ncpp' 
         qui gen double se_`vr'=sqrt(r(Var))
      }

      *Calculating the point estimate as the mean of the effects across draws
      *We don't change the point estimate of the ownership proportions and their differences
      *since we use the actual ones
      foreach vr in `me' `perc' {
         if "`vr'"~="pr_${re}_${as}_`cb'_own" & "`vr'"~="dpus_${re}_${as}_`cb'_tot" & ///
            "`vr'"~="dpde_${re}_${as}_`cb'_tot" & "`vr'"~="pr_${re}_${as}_`cb'_cus" &  ///
            "`vr'"~="dpus_${re}_${as}_`cb'_coeff" & "`vr'"~="dpus_${re}_${as}_`cb'_char" & ///
            "`vr'"~="pr_${re}_${as}_`cb'_cde" & "`vr'"~="dpde_${re}_${as}_`cb'_coeff" & ///
            "`vr'"~="dpde_${re}_${as}_`cb'_char" {
            qui su `vr' in 2/`ncpp', meanonly 
            qui replace `vr' = r(mean) in 1
         } 
      }

      *Saving the dataset with all the probabilities and marginal effects
      qui save "${projf}/Results/ELSA/${as}/boot_alldraweff_${re}_${as}_`cb'.dta", replace

      *Saving the dataset with only the mean magnitudes and their standard errors
      qui keep in 1
      qui save "${projf}/Results/ELSA/${as}/boot_meaneff_${re}_${as}_`cb'.dta", replace
 
      *Saving matrices of marginal effects and percentage variations as a text file

      foreach rs in $prfx  {
         *Number of elements
         local n_`rs': word count ``rs''
         matrix olout=J(`n_`rs'',6,.)

         forvalues k=1/`n_`rs'' {
            local h: word `k' of ``rs''
            matrix olout[`k',1] = `h'[1]
            matrix olout[`k',2] = se_`h'[1]
            matrix olout[`k',3] = olout[`k',1]/olout[`k',2]
            matrix olout[`k',4] = 2*ttail(nobs_`cb' - ${nc},abs(olout[`k',3]))
            matrix olout[`k',5] = olout[`k',1] - invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
            matrix olout[`k',6] = olout[`k',1] + invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
         } 
    
         *Labels of matrix of magnitudes
         *Column labels
         local colmat `rs' se tstat pvalue ci_low ci_high
         matrix colnames olout = `colmat' 
         *Row labels
         matrix rownames olout = ``rs''
         *Displaying and saving the matrix of results
         mat list olout,format(%15.5f)
         mat2txt, mat(olout) saving("${projf}/Results/ELSA/${as}/boot_mefftab_${re}_${as}_`cb'.txt") replace format(%15.5f) ///
            title("Meff for ${as} in ${fnre}, using `cb' quartiles")

         erase "${projf}/Results/ELSA/${as}/temp2.dta"

      }  /* End of loop over magnitudes */        
   } /* End of loop over quartile specifications */
   
   erase "${projf}/Results/ELSA/${as}/temp.dta"         
   }  /*End of loop over countries */    
}  /*End of loop over assets */    

display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"
display "log file for ${as} in ${fnre} closed on $S_DATE at $S_TIME"
cap log close
   
   
  


